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1. Introduction 

Self-organised criticality (SOC) was originally introduced P as an approach to 
understand l//-noise as well as the apparent abundance of power laws in nature, which 
is generally accepted as the sign of scale-invariance. The idea is that under very general 
circumstances driven stochastic processes develop into a scale-invariant state without 
the explicit tuning of parameters, contrary to what one would expect from equilibrium 
critical phenomena jjSj. 

A very large zoo of SOC models has been developed jS], with each model having 
certain special features. However, based on large scale numerical simulations it has 
become increasingly clear that many of the models formerly thought of as representatives 
of entire universality classes or even paradigms for a specific type of model, are either 
not scale- invariant or at least do not follow simple scaling [H El El EI • In fact, models of 
SOC are notorious for slow convergence and deviations from the expected behaviour. 

Some models, however, show all features one would expect from a "self-organised 
critical" system: Consistent exponents and scaling, universality, crossover between 
different classes etc. One of these models is the so-called Oslo model p^, which was 
motivated by an experiment [0]. In a recent paper it has been shown that any 
(small) amount of anisotropy will drive this model eventually (in the thermodynamic 
limit) towards another "fixed point", which is represented by the "totally asymmetric 
Oslo model" (TAOM). This model is solvable directly on the lattice without making 
any scaling assumptions. Consequently, it is not only possible to derive exponents, but 
also to calculate amplitudes of the moments of the relevant observable. 

The TAOM is totally asymmetric in the sense that particles can move in one 
direction only, similar to the totally asymmetric exclusion process (TASEP) |TT| IT^ . 
The TASEP has been solved using a matrix product state ansatz jl2| IT^ . so that it 
seems reasonable to apply similar techniques to the present model. However, there is 
a crucial difference between these two stochastic processes: The relevant observables in 
the TASEP exist on a microscopic timescale, i.e. there is an intrinsic timescale in the 
time-evolution of the microstate of the system. In contrast, in the TAOM the relevant 
observables are obtained by any dynamics which comply with a certain set of rules. In 
that sense, the specific (microscopic) dynamics of the TAOM are irrelevant. This is 
reflected in its theoretical treatment, in that the TASEP is updated homogeneously (all 
sites evolve equally) but the TAOM is perturbed once and is only observed after it is 
fully relaxed (separation of time-scales). 

In the following, the model is defined in terms of rules on a lattice. Using a Markov 
matrix approach it is then solved and exponents and amplitudes derived. After mapping 
it to a reaction-diffusion process as well as various other processes, a more accessible 
continuum theory is described. 
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2. The Model 

The model is defined on a one- dimensional lattice of size L, where each site i = 1,2, ... ,L 
has assigned a slope Zi G {1,2} and a critical slope zf G {1,2}. From a fiat initial 
configuration Zi = 1 and z!^ random, where z![ = 1 is chosen with probability p and 
z"^ = 2 otherwise, the model evolves according to the following rules: 

(i) (Driving) Increase Zi by one unit ("initial kick"). 

(ii) (Toppling) If there is an i where Zi > z"^, decrease Zi by 1 unit, Zi ^ Zi — 1, and 
increase the right nearest neighbour j = i + lhjl, Zj—*Zj + l (charging). A new 
z!^ is chosen at random from {1,2}, where = 1 is chosen with probability p and 

= 2 otherwise. 

(iii) Repeat the second step until Zi < zf ( "stable" ) everywhere. Then proceed with the 
first step. 

During toppling, the right neighbour is charged of course only if it actually exists, i.e. 
j < L, otherwise the toppling site i relaxes without charging another site, so that a 
unit leaves the system. Apart from this boundary condition, the TAOM differs from 
the original Oslo model [Hj in redistributing only a unit to the right, rather than one to 
each side. 

It is important to note that the value of is determined only after a site has 
discharged. Thus, if a stable site i is in state Zi = 1, its value of could be randomly 
chosen in the moment when it is needed, i.e. when the site is charged again. If a stable 
site i is in state Zi = 2, then zf has necessarily the value 2. When this site is charged, 
it will relax to Zj = 2 again and a new random zf is drawn. If that is zf = 1, then the 
site topples again and ends up in state Zj = 1, otherwise it remains in state Zj = 2. 

If all sites are stable, i.e. Zj < zf for all i E [1, L], a configuration is fully described 
by the values of the Zj alone; if Zj = 2 then zf = 2, otherwise zf is random and has not 
yet been used in the dynamics. 

The number of times the second rule is applied, that is the number of topplings, is 
the avalanche size s. The fundamental observable one is interested in is the probability 
density function of these sizes, P{s), which is expected to obey simple scaling above a 
fixed lower cutoff si, 

P{s) = as-^g{s/so) , (1) 

with So = bL^ , Q the universal finite-size scaling function, and metric factors a and b 
[Hj . Equation (jT)) is the definition of the two exponents r and D. It entails that the 
moments (s") of P{s) behave like [TH] 

(s'^) = a(6L^)^+"-"^„ for l + n-r>0 (2) 

with universal amplitudes gn JD] ■ Thus, assuming (0) one can derive r and D from the 
behaviour of any two moments. Below, the exponents 7„ from (s") oc L'^" will be used; 
Equation (j2I) therefore means 

^n = D{l + n-T) . (3) 
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The time series of avalanches, s{t), itself is not Markovian, while the sequence of 
stable configurations of the system, given by the vector {zi, Z2,---, zl), is. Since two 
consecutive stable configurations are not necessarily linked by a unique sequence of 
topplings, the sequence of avalanche sizes is not uniquely determined by the sequence 
of configurations. Nevertheless, in the form of a generating function this ambiguity can 
be built into the Markov matrix operating on the distribution vector of configurations, 
so that the avalanche size distribution can be determined by means of this specially 
prepared Markov matrix. 

2.1. Abelian property 

Put simply, if a model is Abelian jTHj, it means that the order of updates is irrelevant for 
its statistical properties. It is exceptionally simple to see this property here: Firstly, for 
the final state of an individual site there is no difference between a certain number 
of charges arriving at once or arriving sequentially. Secondly, if a site topples, it 
pours particles on its right neighbour, but it will never receive anything back from 
the neighbour. So, if a site at z = 1 has received 3 units, it topples at least twice, but 
for this site it does not make any difference whether it first moves one unit over to the 
right neighbour and waits until all sites to its right have relaxed, or whether it moves 
all units at once, 2 with probability q = 1 — p (namely the probability to have z^ = 2 
after the second toppling) and 3 with probability p. 

In this informal sense, the Abelian property allows the updating to run from left 
to right, completely relaxing each site during a sweep. If there is no toppling on site i, 
the avalanche has stopped and sites j > i do not need to be checked for the toppling 
condition Zj > Zj at all. This procedure makes very efficient Monte Carlo simulations 
possible. Moreover, it defines an activity a^, which is the total number of charges received 
at site i during an avalanche. The activity will be used in Sec. lU 

3. Markov matrix approach 

The tensor product (S> used here is explained in detail in P^l- In particular it has the 
property (provided that A, B, A' and B' have approppriate ranks) 



where stands for the appropriate operator: it is a matrix multiplication if A, B, A' 
and B' are matrices, it is a multiplication of a matrix and a vector if A and B are 
matrices and A' and B' are vectors or vice versa, or it is an inner product if they are all 
vectors. In particular, in the latter case it is 



First, we consider a single site system, which can be in exactly two states, so that its 
distribution of states can be represented by a two-row vector. By convention, the upper 



{A(g)B)Q {A' B') = {AQ A') ®{BQ B') 



(4) 



{a®h){a' ®h') = {aa'){hh') . 



(5) 
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row corresponds to 2; = 1 and the lower row to z = 2. Three matrices are introduced, 
corresponding to the three possible outcomes of a single initial kick. 

The matrix S corresponds to a unit being absorbed, i.e. the site is in state z\ = \ 
and = 2, which occurs with probability q. After the charge, the system is in state 
= 2. Similarly, T corresponds to a single toppling due to the charge and U corresponds 
to a double toppling: 

-(^:) -co 

In the following, the aim is to find an expression for the moment generating function of 
the avalanche size distribution. To this end, each matrix is multiplied by an appropriate 
power of X, so that evaluating at x = 1 gives the usual Markov matrix of this process, 
and deriving by x before evaluating at x = 1 multiplies each process by the number of 
topplings occurring in it, and similarly for higher order moments [T7] . 

It will be motivated only a "posteriori that a dissipation process is required, say with 
probability < e < 1; this process corresponds to charging without changing the state, 
i.e. an identity operation 1, the latter being the 2x2 identity matrix. The resulting 
single site operator is therefore 

with 5 = 1 — e. The eigenvectors and eigenvalues of this matrix are found to be 

(eA(x)l =(i,l) 

(e.(x)l =(-f,p) 

where Oi acts on bra-vectors (| from the right and on ket- vectors |) from the left. The 
vectors are normalised such that 

(ealCfe) = 5a,b (9) 

with ba), denoting the Kronecker delta-function. In order to distinguish vectors of 
different size, in the following they are often marked with an index L to indicate a 
size 2^. 

Ol(x) is the operator which adds a unit on site i = \ and relaxes the entire lattice 
of size L. It is a matrix of size 1^ x 2^ and defined as 

Ol(x) = el®^ ^b(S® + xT ® Ol_i(x) + x^U ® 0\_^[x)) (10) 

again with a dissipation rate e, leaving the state unchanged. The bracket multiplied by 
5 consists of three terms: The first term charges the site without toppling and leaves 
the rest of the system unchanged by operating with the identity The second 

term corresponds to a single toppling, which charges the remaining system of size L — 1 
once. This term is derived using the identity 

(t ® 1^(^-1)) (1 ® Ol-i{x)) =T® Ol-i{x) . (11) 



xp 

1 . 

—X 

1 



A = e + x5 
/i = e 



(8) 
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The third term is a double topphng of the site, giving rise to a double charge of the 
remaining system. 

The Abelian property mentioned above (Sec. I2.1|l can be expressed as the 
commutator for two charges on a system of size L, one at site i = 1, the other one 
at site 1 + L — L' with L' being the size of the subsystem starting from the site receiving 
the second charge, 

Ol{x) ® Ol,{x)) = Oz.'(x)) Ol{x) (12) 

where of course L > L'. The tensor multiplication used on O^/ and also in (fTT|l ensures 
that both matrices have the same rank; they are "filled with identity" where they do 
not act. Equation (fT2|l simply states that it does not matter for the statistics whether 
the leftmost site of a right subsystem of size L' in a system of size L is charged first, 
followed by the leftmost site of the entire system, or vice versa. Due to the asymmetry 
in the dynamics, it is clear that a system of size L, initially charged at site i, has the 
same statistics as a system of size L — i + 1, charged at its leftmost site. It might be 
interesting, however, to formally prove Equation (fT^ . which should be easily feasible 
using established methods [THl IT3]. 

The distribution of states at time t is the vector \Pt)i, which has rank 2^, each row 
corresponding to the probabihty for the system to be in the state encoded by that row. 
The encoding follows from the row ordering convention introduced above and the use 
of the tensor product in (fTUjl . 

For X = 1 the operator Ol{x) is simply the Markov matrix acting on |-Pt)j^, 
producing the distribution of states at time t + 1 [Hj 

|P,+i)^ = 0^(1) \Pt)^ . (13) 

There exists at least one eigenvector with eigenvalue 1, which is therefore a stationary 
distribution. If the eigenvectors represent a complete basis and the modulus of all other 
eigenvalues is less than unity, this stationary distribution is unique and reached by 
any initial distribution. The stationary distribution, denoted |0)^, is the focus of the 
following calculations. It is shown below that it is unique. 

One very important bra-eigenvector with eigenvalue 1 of Ol{x = 1) is 

(0|^ = (M^_^) (14) 

2^ times 

by normalisation. As has been indicated above, for general x, the operator Ol{x) 
becomes a moment generating function of the avalanche size, if sandwiched between 
{0\j^ and the stationary distribution: 

QL,ni^;e)^{0\,Olix)\0)^ (15) 

This can be seen from (fTn|) containing an x for every toppling. When the operator acts on 
a distribution, for each transition from one state to another a power of x corresponding 
to the number of topplings enters and is multiplied by the probability to be in the 
initial state (given by the initial distribution) and the transition probability given by 
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the transition matrix. The function QL,„(x;e) for general n is then the generating 
function for avalanches caused by n = 1, 2, . . . initial kicks. In particular 



X 



dx 



QL,l{x\t) . 



(16) 



x=l 



The aim of the following calculations is to find the generating function Qi^i{x\e) 
or at least the moments generated by it. 



3.1. General eigenvectors and eigenvalues ofOL{x) 

It would be very helpful if 0^(0;) could be written in the form 

OL{x) = 'Y:m)L^L,m\L ^ (17) 

i=0 

where {i{x)\]^ denote the left hand and \i{x))j^ the right hand eigenvectors of Ol{x) with 
eigenvalues Xi^iix) and i = ... 2^ — 1. A priori it is not clear whether these vectors 
actually exist. In the following they are constructed and it is shown that setting e = 
leads to fundamental problems. 

Assuming that \i)^_i is an eigenvector with eigenvalue XL-i,i of Ol_i(x), the 
definition of Ol(x), Equation (fTUj). gives for an arbitrary vector |e)^ 



Oi(s) (|e)i ® |z)^„i) = [{el + 5[S + xTXL-l,^ + x''UXl_,^,)} \e)^] ® \i) 



'L-1 



where |e)-,^ contains two elements such that \e)^ ® K)^,-! is a vector of 2^ elements. The 
matrix in the curly brackets is simply Oi(xAL_i,i). So, if \e)^ is either \e\{xXL~-i,i)) or 
|e^(a;Ai_i_j)) from Q, then |e)^ ® is an eigenvector of Ol{x) with eigenvalues 

e + 6{xXL-i,i) or e. Thus, based on Q, one can write the eigenvectors of Ol{x) as 



K)l 
(^II 



i + 2^-1 
i + 2^-1 



\ex{xXL~i,i)) ® \i)L-i 
{ex{xXL^i^i)\ O {i\^_^ 
\ef,ixXL^i,i)) O \i)L_i 



(19) 



and the eij 


2;envalues 


as 






AL,i = e 




^L,i+ 


2L-1 = e 


both with 


1 = 0,1,. 


..,2^-1- 




|0)i 


= \(^x{x)) 




(Ok 


= (eA(a;)| 




ll)l 


= \ef.{x)) 




(111 


= (e/.(a;)| 


and the ei; 


2;envalues 


as 




Ai,o 


= e + x6 




Ai,i 


= e 



(20) 



1. To start the hierarchy, one defines 



(21) 



(22) 
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Now it is clear why the quantity e was necessary: For e = all but one eigenvalues 
vanish, which can be seen from the hierarchy of eigenvalues obtained by iterating (pUj) . 
Therefore, if the vanishing eigenvalue of L — 1 is plugged into (caI or (e^| according to 
(fT^. the result is undefined, as can be seen from Q, so that the bra-eigenvectors cease 
to exist. 

The fact that all but one eigenvalues vanish for e = is very deceptive. Assuming 
that any initial condition \P) can be written in terms of the eigenvectors of 0^(1), say 
X^Oi 1^), this suggests 0^(1) \P) = |0). This, however, is wrong, because for vanishing 
e the operator Ol{x) cannot be written in the form ()17|) for L > 1. And it must be 
wrong, because, for example, kicking an empty system once cannot make it produce the 
stationary distribution. 

If the eigenvectors of Ol-i are linearly independent, then one can show, using the 
construction (fTn|l . that the eigenvectors of are linearly independent as well, provided 
that \ex{xXL-i,i)) and \e^{xXL_i^i)) are linearly independent. This is not the case for 
e = (see the ket vectors in ^ with x = 0) and this is the basic reason why e 7^ 
is needed for the time being. However, for any e 7^ one can apparently construct a 
diagonalising matrix for O^. Thus, it can be written in the form (|T7jl . Especially, the 
eigenvectors have the property (by induction) 

{tix)\jix))^ = S,, (23) 

and as all 2^ eigenvectors are linearly independent, they must span the whole space so 
that 

K)l (^II = 1 ■ (24) 

i=0 

In the form ()17j) the operator can now be applied to a stationary distribution to give 
QLA^;e) = '^' mx)),Xl,{tixm^ (25) 



i=0 



3.1.1. The stationary distribution From (|T!Hl the stationary distribution can be derived 
immediately. It is the eigenvector with eigenvalue 1 of 0^,(1). Setting x = 1 in ()2Up 
it is clear that Xi^i = 1 requires a Xi-ij = 1, which, together with (j22I), gives the 
unique A^^o = 1 provided that e < 1. If e = 1, then all eigenvalues are 1, but still 
all eigenvectors are linearly independent and therefore span the entire space, so that 
all initial distributions are stationary. This is not surprising because e = 1 simply 
means that any added particle immediately dissipates from the system, so that adding 
a particle is in fact just the identity operation. 

For < e < 1 the stationary distribution is unique and all other eigenvalues have 
modulus less than 1. The eigenvector corresponding to eigenvalue Al^o = 1 is, according 
to (Uni), 

(0|, = (e,(l)r^ = (l,l)«^ (26a) 
= |eA(l))^^ = M (26b) 
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which is consistent with the notation for the stationary distribution and the 
normahsation eigenvector introduced in ()15|) and (fT^ . The last hue, Equation ()26b|) . 
indicates that the stationary state is a product measure, i.e. a state at one site does not 
depend on the state on any other site. In fact the spatial correlation function of sites 
{ii,i2, ■ ■ ■} can easily be calculated by "dressing" the states of the sites by appropriate 
powers of a variable Xi, in order to obtain the generating function of the correlators. 
The function 

Cix^,x„ ...,x,) = {0\J ]( ^""A ] = liiP^^ + (27) 

\ J \ qx2 J V I^L J i 

is the generating function of the state-correlators, where state 1 stands for 2; = 1 and 

state —1 for z = 2. The states have the useful property that the joint contribution of 

two sites is 1 if both sites are in the same state and —1 otherwise. The average state is 

obtained from 

d 



ln(C(xi,X2, . . .,xl)) = 

Xl,...,XL = l * 



C{xi,X2,...,Xl) = p-q (28) 

Xl,...,XL=l 

Correspondingly, the connected two point correlation function of sites i and j is given 
by 

for i = j (29a) 



d d 




otherwise (29b) 



In (C(xi,X2, . . ■,xl)) = 

Xl,...,XL = l 

This confirms the absence of correlations and is fully consistent with the expected 
variance of the state. 



3.2. The hierarchy of generating functions 

Using (j2Sl), one can now calculate the generating function QLn{x;e), by plugging the 
hierarchy of eigenvectors (fT^ and eigenvalues ()2U|) into ()25|) and using the properties of 
0, see Equation For n = 1 it is 

2^-1-1 

g^,i(x; e) = ^ xS\l-i, m^., {0\exiXL-i,)), (^|0)^_i (eA(Ai_i,)|0)^ 

i=0 

+ '^'e(0|z(a;))^(z(x)|0)^ 

i=0 

where the term proportional to e comes from the e in every A^^j (see Equation ()20|) ). 
From (j24j) it is clear that the last sum gives 1. The two projections give 



P 

(0|eA(AL_i,i))^ = \L-i,iP + q and (eA(AL_i,i)|0)^ = — h q 



so that 
Ql,i{x; e) 



(30a) 



e + ^ x6XL-i,i (0|i)^_, (2|0)^_, (^p^ + q^+pq (^x\L-l,^ + ^^^^^^ ) 
e + (5((]9^ + g^)x(5L_i,i(x; e) + pq{x^QL-i^2{x\ e) + QL-ifl{x\ e))) (30b) 
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where Equation (^3]) has been used in the last hne. Of course, the generating function 
QL,n{x', e) is defined, p5|) . for all e and therefore one can take the limit e 0. This limit 
should not cause any problems, as e has only been used to construct the eigenvectors. 
In fact, the limit must be identical to setting e = in (jHUj) . as can be shown from ()3Up 
by induction in L. This finally gives 

Ql,i(x; 0) = {p^ + q^)xQL-i,i{x; 0) + pq [x^Ql-iA^] 0) + Ql^i,o{x; 0)) , (31) 

where QL-ifi{x; e) = 1 by Equation (fTH|) . consistent with Equation (j2Sl)- In fact, the 
calculation above can be generalised: 

i=0 j=l \J J 

X (0|z)^_, (0|eA(AL_i,))i (^|0)i_i {ex{XL-imi (32) 
+ '^'e"(OK(x))J^(x)|0)^ 
Again, all sums can be written in terms of QL-i,nix; e) plus e": 



X 



{{P^ + q^)x^QL~i,j{x; e) + pqx^^^QL-i,j+i{x] e) + pqx^ ^QL-i,j-i{x; e)) (33) 
For vanishing dissipation this simplifies to the central result 

Ql+i,„(x; 0) = x" (dQlAx; 0) + D{xQL,n+i{x] 0) + x-^QL,„_i(a;; 0))) (34) 

with D = pq and D = + = 1 — 2D. Equation ()34j) is closely related to a diffusion 
equation. The boundary conditions are QL,o{x;e) = 1 for L > 1 as mentioned above 
and QL=o,n{x] e) = 1. The latter comes from a direct evaluation of for L = 1, which 
is identical to for QL=o,n{x; e) = 1. 

3.3. Solving QL,n 

There is no general solution for ()34j) known to the author. However, one can solve it 
order by order in derivatives by x at x = 1, i.e. calculate every individual moment, see 
(fTB|) . In the following, the notation 

QL+l,n = QL+l,n(l; 0) (35a) 

Ql+i,„(x;0), (35b) 

x=l 

etc. is used. One finds for n > 1 



^^+1'" dx 



QL+l,n = DQi^n + D{QL,n+l + Ql.^-i) (36) 

which is solved with the boundary conditions introduced above by QL,n = 1- Of course, 
this is just normalisation. Using this result the next derivative is 

QL+l,n = n + DQ'L,r. + ^(^L.n+l + Q'l^u-i) (37) 
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with the boundary conditions Q^_,_i„=o ~ Q'L=on — 0- The solution of (jHTjl can 
easily be guessed as 

Q'L,n = nL . (38) 

This is not surprising, because it says that the average number of topplings occurring 
in the system per n kicks is riL. That is obviously true, because every unit added must 
leave the system by travelling through the entire lattice. 

The next order is the first non-trivial one. The difference equation then reads 

Q'Ui^n = {n' + 2D){2L + l)-n + DQl, + D{Ql^^, + Ql^_,) . (39) 

Introducing 

Q'Ul,n = SL+l,n + QU,n (40) 

with 

L-l 

SL,n = J2n^-n + 2in^ = -nL + n'^L'^ (41) 

which has the useful property SL+i,n — Si^n = —n + 2n'^L + and SL,n+i + Si^n-i = 
2SL,n + 2L^, one finds after some algebra 

' ' l\ ( I \q\ 



Q'i, = Y.mL-l?j:{%''-'"'q'-' 



(42) 



^ \m J \m + IJ p J 
In order to analyse the asymptotic behaviour for L ^ oo, one writes 

Qli = 2^(^-0 V(0 (43) 

1=0 

with 



y \m + IJ p J 

The binomials, which can be approximated by a Gaussian, are treated identically and 
the summation can be written as an integral, so that finally 

Q'U - I''' 2D(L - If-^ (Si^q) + S{^)^ (45) 
where £{x) = 2 f^ dzex^ii^—z^) / ^ . In leading order, this turns out to be 

which is according to (j40|) also the leading order of Q'^+i „ and therefore the leading 
order of (s^), see (llfjj) with This is perfectly confirmed by numerical simulations 

of the model. 

The two exponents 71 = 1 (see Equation (fHHjl ) and 72 = 5/2 pHjl lead together 
with © to r = 4/3 and D = 3/2. 
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4. Reaction-Diffusion mapping 

It is possible to map the model onto a very simple reaction-diffusion process of the form 
A + A ^ A pni- To this end, the configuration of the lattice is described by the thick 
line shown in Figure Q The line consists of segments, which can either point up or down 
by an angle of 45°. If the line corresponding to the ith site goes up, it indicates that the 
ith site is in state z = 1, otherwise the line goes down indicating the state of the site to 
he z = 2. According to (j26bjl the configuration of the lattice (in the stationary state) 
after an avalanche is a product state, where a site is in state z = 1 with probability p 
and in state z = 2 with probability q. Thus, the thick line is in fact the trajectory of a 
random walker with drift corresponding to the difference p — q. 
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Figure 1. The thick, full line shows the configuration of the lattice after an avalanche 
has passed through. Each up or down-pointing segment corresponds to a single site, 
the position label of which is shown under the dotted line. A segment pointing upwards 
corresponds to a site being in state z — 1 (with probability p, see Equation Ij26b|l ). a 
segment pointing downwards corresponds to state z = 2 (probability q), as indicated. 
The dashed line corresponds to a "toppling trajectory" as explained in the text. 

The avalanche itself, on the other hand, is a random walk with the same 
probabilities. One can see that by considering the activity a^, which is the number 
charges received at site i during an update-sweep as described in Sec. 12.11 From 
site to site, the activity can either remain constant or change by 1 up or down. 
Apparently, ai = 1 is the driving. If a site receives Oj charges and changes state by 
/S.Zi = Ziit) — Ziit + 1) , then its right neighbour receives Oj+i = + A^j charges. If 
Azi = 0, then the vertical distance between two consecutive configuration trajectories 
(as shown as thick and dashed lines in Figure does not change. If, however, the new 
configuration of site i has an increased value Zi{t + 1) > Zi{t), the activity goes down, 
because Azi < 0. The only way to increase Zi is to go from state 1 to state 2, i.e. the hne 
segment of the former configuration points up (probability p), while the line segment 
of the new configuration points down (probability q), so that the gap between the two 
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trajectory decreases, see Figure [T] at the dotted line 3. Similarly, if the activity goes up, 
then Azi > and the gap increases. 

After the activity vanishes, the profile of the new configuration remains unchanged 
compared to the former, i.e. the gap between the two configurations is a constant. In 
fact, if the gap was initially 1 and goes up and down by 1 as described above, then 
the gap will be as soon as the activity vanishes. This is exactly what is shown in 
Figure^ The thick dashed line shows the new configuration and its distance to the old 
configuration is the activity during the avalanche. This avalanche occurs within the 
configuration shown as a thick line, initiated by a single kick. Initially the gap is ai = 1. 
If the dashed line would go down immediately on site 1, the site would have "absorbed" 
the initial unit and would be in state z = 2 (i.e. a segment pointing down). Instead, in 
the example, it goes up twice; first just like in the old configuration so that the activity 
does not increase, and then in the opposite direction to the old configuration so that 
the activity increases by 1. On site 3 and 4 it goes down twice; the toppling on site 
4 is particularly interesting. Here, initially the activity is 1, i.e. the site has received 
one unit. But the site is in state z = 1, so it absorbs the unit with probability q, 
corresponding to the probability of the dashed line segment to point downwards. 

The activity is measured half a unit left of each site as the distance between old 
and new trajectories, which, in turn, is measured in such units that the vertical distance 
between two circles (in Figure is 1. The reason for the shift is that one wants to 
measure how many charges have arrived at a site, not affected by the value of the 
resulting activity. 

To repeat this important point, the trajectory of an avalanche becomes the 
configuration for the next avalanche, i.e. the thick dashed line in Figure C] becomes 
the thick solid line for the next avalanche. 

One can calculate the probability of the changes of activity explicitly: The new 
segment goes up with probability p and down with probability q, the same applies to 
the old segment. Thus, they point in the same direction (no change of activity) with 
probability + g^, the gap widens with probability pq and shrinks with qp. Hence, the 
gap between the two trajectories is in fact a symmetric random walk, even though the 
individual trajectories might have a bias, according to p — q. 

As described above (see Sec. EI), the avalanche size is measured as the number of 
topplings. For convenience, one can define it as the number of charges, which makes 
hardly any difference, because the number of topplings of site i is identical to the number 
of charges on site i + 1, unless i = L, simply because there is no site L + 1; similarly for 
L = l. 

In the following, we will consider the number of charges as the avalanche size, 
because the total number of charges is simply the area between two of those trajectories 
described above, namely the sum over all activities. From this it is also clear that the 
avalanche size is actually uniquely determined by the initial and the final configuration, 
with initial activity ai = 1. 
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4.1. Relation to other models 

Before the above identification of the process as a random walker is cast into an 
continuum problem and subsequently solved, it is worth pointing out other models 
which are closely linked to the present one. 

4-1.1. Anisotropic BTW model Dhar and Ramaswamy jJOj developed an anisotropic 
variant of the well-known BTW sandpile model PP , which is now known as the directed 
sandpile model. This model, however, is situated on a 1 + 1-dimensional lattice and the 
annihilating random walkers represent the contours of the compact area covered by an 
avalanche. The randomness here comes solely from the randomness of whether a site 
charged by particles from toppling sites topples in turn. An equivalence to a variant of 
directed percolation has already been pointed out in see also 

Kloster, Maslov and Tang have studied a stochastic directed sandpile model, 
which was originally proposed by Pastor-Satorras and Vespignani This model is 
closely related to the one presented in this paper, even though it is also situated on a 1+1- 
dimensional lattice. The authors find the same exponents by scaling arguments. The 
mapping to the two-dimensional react ion- diffusion process presented above, questions 
their assertion that their model is in a different universality class than the model by 
Dhar and Ramaswamy. 

For these models it is fairly obvious how to extend them systematically to higher 
dimensions. Using scaling arguments in conjunction with some simplifying assumptions, 
Paczuski and Bassler arrive at a general expression for the value of the exponents of 
this model in higher dimensions. Unfortunately, it is not so clear how to generalise the 
model studied in this paper to higher dimensions, because it is unclear how to generalise 
the driving and what boundary conditions to apply. 

5. Continuum solution 




t outflow 



Figure 2. The area under the trajectory (hatched) is the avalanche size. The two 
filled circles mark the starting point (0, xq) and the end point (t, x) 

Having mentioned already the mapping to an annihilating random walk, the 
continuum description is straight forward. To this end, the quantity lijnit^x^xo) is 
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introduced. It quantifies the properties of a random walker along an absorbing wall. 
For n = it is the probability density of random walkers at time t and height x over 
the absorbing wall, starting at height xq, which is xq = 1 for a single initial kick. Here, 
t takes on the role of the horizontal (continuous) position between t = and t = L in a 
picture like Figure To motivate the following calculation, one imagines a large set of 
trajectories of random walkers along the absorbing wall from t = 0, x = Xq to t and x. 
The set of areas under the trajectories, as exemplified in Figure |21 is then {si{t, x; Xq)}, 
where i is indexing the elements in the set. ({s"'(t, x; xq)}) is the average of the nth 
moment over this set. Now one can express the time-evolution of this average as the 
sum of three contributions of the three processes of up, down or straight movement of 
the random walker. Thus, up to terms of order At Ax (see caption of Figure Ej) 



where each term corresponds to a process like the one shown in Figure El The 
multiplication by iljQ{t,x;xo) is necessary in order to weight each of the ensembles for 
each contribution properly. For example, there might a much larger contribution from 
below, even though on average the moment at this position is smaller than at the other 
positions. 



Figure 3. A new segment (hatched area) is added to the currently considered path, 
increasing all areas in the ensemble {si{t, x—Ax; Xq)} by xAt+0{AtAx) and producing 
a new ensemble {si{t + At,x;xQ)}. The example shown corresponds to (|47c|) . which 
starts at (t, x — Ax). Starting points of other contributions are shown as empty circles. 
The coordinates of the two black points are given in the form {t, x). 



ipQit + At, x; Xq) ({s"(t + At, x; Xq)}) 

= pqipoit, X + Ax; Xo)({(s(t, x + Ax; Xq) + xAt)"}) 
+ (p2 + q^)^Q{t, x; xo) ({(s(t, x; Xq) + xAt)"}) 
+pq'iljo{t,x - Ax; Xo)({(s(t, x - Ax; Xq) + xAt)"}) 



(47b) 
(47c) 





At 



Defining 




(48) 



one finds in the continuum limit of (P7j) (keeping DAt/Ax"^ constant) 

dt1pn(t,X] Xo) = DdliJn{t,X]Xo) + Xmpn-lit , X] Xq) 



(49) 
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where D = pq again|. The boundary conditions for n = are observed immediately 
and transfered to ipn using (jlH)) by noting that {{s(t,x; Xq)"}) is non-divergent, so 

limipnit, x; Xq) = 6nfiS{x — Xq) (50a) 

Mt, 0; Xo) = (50b) 

and the PDE (jlHl) is to be solved for x G [0, oo[. 

The avalanche sizes are measured from avalanche trajectories which have died out 
or reached the end of the system. Thus, the averages measured in the model are taken 
from the random walkers which have reached the absorbing wall or did not do so until 
a cutoff time t. Therefore the nth moment observed is 

rt roo 

(s")(t;xo)=/ dt'in{t']Xa) + dx'ipnit^x'-.xo) (51) 



where the first integral runs over the "outflow", j„(t, x = 0; Xq) = —Ddx\x=o'ipn(t, x; Xq) 
and the second over the contributions at cutoff time (see marks in Figure |2|). (<s")(t; xq) 
denotes the nth moment of the avalanche size (measured as the number of charges) for 
a system of size t starting with xq initial charges. Using one has 

rt roo 

(s")(t;xo)= / dt' dx' x' nip n-i{t' ■, x' ] Xo) (52) 
JO Jo 

The dimensionless form of ip is given by 

1 [ x'^ \ " 

ipn{x,t]Xo) = — \^] i'niy,r) (53) 
Xo \D J 

with y = x/xq and t = t/ (xl/D). The propagator G{y, r; yo) is easily obtained from a 
mirror-charge trick. 



1 



y-VQ v+yp 



G{y, t; yo) = -j= ( e 4. - e 4r ) , (54) 



and i'o{y,T) = G{y,T; 1), i.e. 



^o(y,^) = -^=<^ sinh ) . (55) 



One might be inclined to transfer the problem into fc-space, which, however, does not 
simplify the problem because of the boundary condition ()50b|l . The expression 



i^n{y,T)= dr' dy'ny'ilJn-i{y',r')G{y,T - T';y') (56) 

JO JO 

is the formal solution. Rescaling the arguments of ipn by powers of /i one finds 

^n{^/liy, jJ-r) = ^^1'^ / dr' \ dy'ny'ipn~i{y/J^y', I^T')G{y, r - r'; y') . (57) 



From il)n-i{y/l^y,^-'r)= /i""-i^/'„^i(?/, r)^then follows ipn{y/J^y, fJ'T) = /i""-i+^/2^„(2/, r). 
Thus, starting with Tpo^^/Jly , ^it) = ^°'°ipo{y,T) one has apparently 

MVJ^y, t^r) = /it"+"V„(y, r) . (58) 
I It is interesting to note that this can be written using a generating function 'I'(t, a;; a;oi with 



= x^^f + so that indeed ^ "ifit, x; Xq, £,) ^ tpn{t, x; Xq) 
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MVJ^y^l^r) = -^e-'^(l + -(^-^\+...\ . (59) 



Unfortunately the scaling behaviour of ipQ is a bit more complicated. Nevertheless, it 
can be expanded for large fi, or actually large fir, as 

1 I ^-y^ ( y 1 ( y_' 

V2r /i V48r3 Sr^ 

Bearing in mind the necessity of large fir one can now apply the scaling argument ()58|] 
order by order in fi since Equation P9|l and its dimensionless counterpart are linear. 
From (jnni) it is ao = — 1 for the leading order, ao = —2 for the first sub- leading order 
and so on. 

Equation (jKH|) immediately translates to {s^) using (jK^ and to leading order 
one finds 

(S") {fit; Xo) = ^(3/2)n+l/2+-0 ^,n^ ^f, Xo) + . . . (60) 

Assuming (^, from (j2I) with t taking the role of L it follows that D = 3/2 and 
D{1 — r) = 1/2 + oq, i.e. for = —1 one has r = 4/3. The next order correction is 
D = 3/2 and r' = 2. 

Of course, it is also possible to calculate the leading orders of {s^) exactly. Because 
of (j60|l . one needs to calculate (s")(/it;xo) for one value of t only. The simplest choice 
is to set if: = x^/D, which gives {s){xq/ D]Xq) = x^/ D for n = 1, i.e. 

{s){t;xo) = Xot (61) 

which is exactly (j38|l {n in (j38|l corresponds to xq here and L in (j38|l to t). This is 
actually surprising, because (jHT|) is only the leading order and corrections are expected 
from higher orders. However, it turns out that in fact all higher order corrections cancel. 
Indeed, remarkably 

dxxe~^ (^j - 2e"3t sinh (^^^^ = . (62) 

even though x/t is only the leading order of 2exp(— l/(4t)) sinh(a;/(2t)). Especially 

' , f x^ X \ . . 

dxxe-^^ (i^-^j=0- (63) 

According to the next moment is 

{s'){f,t;xo) = 2f,'/' drj^ dy^yMy.r) (64) 

the leading order of which can be determined using the leading order of ipi, 

r , , 1 v'^ v'^ 1 / y-y'. y+y' . \ 

tlj^(y^T)= / dr' dy'—=e~^ , . e W^) - e + . . . (65) 

"'o Jo \/t't[ ^'T JA-Kij — t') ^ ^ 

which gives the leading order of (s^) 



32 



it- xo) = ^t'/'^Dxl + 0{t'") (66) 

identical to P6|l . Higher orders become very tedious, so that numerical evaluation seems 
to offer the better option. 
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6. Discussion and Conclusion 

The results above represent some of the few exact result for sandpile-like models: 
Equation ()38|) and Equation (jlUj) are the exact leading orders of the first two moments 
of the avalanche size distribution without making any assumptions about scaling 
behaviour. The conclusion that r = 4/3 and D = 3/2 can only be drawn by either 
assuming (H)), or by accepting the continuum result and using the uniqueness of the 
distribution inferred from its moments. § 

The method introduced in Sec. El is not restricted to sandpile-like models. The 
underlying idea is to use a Markov matrix not only to evolve the state distribution, but 
also to calculate the moment generating function of the relevant observable. In order to 
obtain the finite-size scaling behaviour, its set of eigenvectors is generated recursively. 
From this recursion relation one can then develop a (discrete) PDE like (jH^), which can 
subsequently be used as a starting point for other techniques. In a two-dimensional 
variant of the present model, this recursion relation is much more complicated to obtain 
and might require the use of a matrix product state ansatz ^Hj- Nevertheless, it seems 
promising to apply the approach to more complicated processes, such as the TASEP 
and recent variants pS], for which there is no solution known yet. 
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